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Monte Carlo simulations have been used to study the phase diagrams for square Ising-lattice gas 
models with two-body and three-body interactions for values of interaction parameters in a range 
that has not been previously considered. We find unexpected qualitative differences as compared 
with predictions made on general grounds. 
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Experimental studies of phase transitions in ad- 
sorbed monolayers have resulted in examination of order- 
disorder transitions in lattice gas Ising models which rep- 
resent the occupation of the periodic minima in the sub- 
strate potentialfll, '^l. Such models, usually containing 
two or more competing two-body interactions, have been 
studied by Monte Carlo simulations [3, 0] which have de- 
termined the location and nature of the resultant phase 
boundaries. Typical ordered phases which are found for 
square lattice models with near-neighbor coupling are 
shown in Fig. [T] along with low density and high den- 
sity disordered states termed lattice gas (L.G.) and lat- 
tice liquid (L.L.), respectively. Experimentally observed 
asymmetries in phase boundaries as a function of cover- 
age can be explained by lattice gas models if three-body 
interactions are introduced as well. Theoretical studies of 
adatom-adatom interactions find surface mediated many- 
body couplings 0. A square lattice model with first- and 
second-neighbor two-body coupling and weak three-body 
interactions was investigated [y] in an attempt to clarify 
the phase transitions for H on Pd(OOl) and predictions 
were made for the case of even larger three-body cou- 
pling. Monte Carlo simulations showed that the inclusion 
of three-body interactions did make the transition asym- 
metric about 50 percent coverage and could even force 
the tricritical point on one side of the phase boundary 
to zero temperature. Such an asymmetry can also be in- 
terpreted in another form of non-additive interactions f^, 
and a recent Monte Carlo study shows the effect of vari- 
ous cases of strength of interactions Q. 

In this paper, we present the results of an investigation 
of the model proposed in Ref fo*! with moderate to large 
three-body interactions using Monte Carlo simulations. 
In the next section, we shall review some appropriate 
background and in Sec. Ill we present our results for 
two different values of interaction parameters which yield 
qualitatively different phase diagrams. We conclude in 
Sec. IV. 



A lattice gas model is a collection of atoms whose po- 
sitions may take on only discrete positions which form 
a periodic array, in this case a simple square lattice. A 
configuration of this lattice is defined by site occupation 
variables Ci where ci = 1 if site i is occupied and q = 
if the site is empty. The Hamiltonian which we use in- 
cludes interaction (/3„„ between nearest-neighbors, (pnnn 
between next-nearest-neighbors, and ipt between neigh- 
bors on a triangle inscribed within a square made up of 
nearest- neighbors: 

U - flNa = - (e + ^) ^ Q - ipnn ^ CiCj 
i i^j 

- Vnnn'^CiCk - ^Pt ^ QCjCfc, (1) 

and the coverage of the lattice is given by 



(2) 



This model may be transcribed to the Ising model by the 
transformation to spin variables ct = 1 — 2ci, thus giving 
rise to the Hamiltonian 
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H'^^Gi — Jnn ^ CrjfTj 



Jt <J^aJak, (3) 

i^k i^J^fc 

where the Ising model interaction parameters are related 
to the lattice gas couplings by 
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Jnn — -^^nn + -^ft 
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FIG. 1: (a) Schematic view of 
the (100) surface of a substrate 
whose periodic potential pro- 
vides a square lattice of pre- 
ferred adsorption sites. The in- 
teractions used in this study are 
shown schematically as straight 
lines between adatoms which 
are represented by the filled cir- 
cles, (b) Unit cells of the or- 
dered overlayer structures dis- 
cussed in the text. 



Jt = -^Vt (6) 



H ^ + ^l) ~ ipnn - '^nnn-'^'^t (7) 

The magnetization is then related trivially to the cover- 
age 

m = 1 - 26*. (8) 

Because of technical considerations, it is generally easier 
to carry out simulations in the magnetic (Ising) represen- 
tation, and additional symmetries often become obvious 
in this approach. For example, in the Ising representation 
it is easy to see that the phase diagram must be symmet- 
ric in the absence of three-body interactions. Throughout 
the remainder of this paper we shall normalize all quan- 
tities by the nearest-neighbor coupling J„„ and define 
R = Jnnn/Jnn and Rt = Jt/ Jnn- As a function of the 
relative strength of the three-body interactions. Binder 
and Landau suggested four different schematic phase di- 
agrams (shown in Fig. [2|) to describe the range of possi- 
ble behavior due to the inclusion of three-body coupling, 
and using Monte Carlo simulations verified the nature of 
the phase diagrams with weak three-spin coupling shown 
in the upper portion of Fig. [21 The remaining two di- 
agrams were not based on any explicit calculation but 
were guessed as extensions of what was then believed to 
be the field-dependent behavior in the absence of three- 
spin interactions. Our views of the correct behavior in 
this latter case have changed, and given the complex- 
ity of phase diagrams in other two-dimensional systems 
with competing interactions, the predictions should be 
regarded with care. 

At T = the energies of each ordered state as well as 
the lattice gas and lattice liquid phases may be calculated 
without difficulty and the intersections of the energy vs. 
chemical potential lines locate the transitions between 
neighboring phases. Of course, this is also valid in the 
Ising representation, and in Fig. [3] we show the ground 
state phase diagrams as a function of Rt and H which 
we obtain for two different values of R which are in the 



parameter region discussed by Binder and LandauQ but 
not actually simulated. 

One interesting feature of this diagram is that the 
(2x2) state is actually degenerate in that either alter- 
nate rows or alternate columns may be shifted randornly 
by one lattice constant without any cost in energyjl, 
This "row-shifted (2 x 2)" state has been seen before and 
the nature of the finite temperature transition to the dis- 
ordered state is a matter of some disagreement. We be- 
lieve that as long as the ground states do not change, 
the specific choice of parameters is not important and 
we have simply chosen values for which the c(2 x 2) and 
row-shifted (2 x 2) phases are stable over relatively large 
ranges of fields in Fig.[3]at which multiple phases become 
degenerate. 

We have used the parallel tempering algorithm with 
GPU accelerated techniques to study the thermody- 
namic properties of this model from which phase dia- 
grams can be deduced. Spins were placed on L x L 
square lattices with periodic boundary conditions and 
were flipped using a Metropolis transition probability. 
Typically, 10^ to 10^ Monte Carlo steps are used to col- 
lect data for each run and 3 to 6 independent runs are 
taken to calculate standard statistical error bars, and 
in all the plots of data and analysis shown in follow- 
ing sections, if error bars are not shown they are always 
smaller than the size of the symbols. Lattice sizes from 
L = 32 to L = 300 were simulated and the data were 
interpreted within the context of finite size scaling [lo|. 
Most of the simulations were carried out on a GeForce 
GTX285 graphics unit. In addition to internal energy, 
specific heat, and magnetization, we calculated order pa- 
rameters, e.g., ?Tic(2x2)i '712x1, etc., for the various or- 
dered states shown in Fig. [U and 4th order cumulant U 
is defined in terms of the order parameter accordingly. 



III. RESULTS 

A. i?= 1/4, i?t = -1/4 

Bulk properties such as the specific heat peak, temper- 
ature dependence of the 4th order cumulant of the order 
parameter, etc., were used to determine the location of 
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FIG. 2: Possible schematic temperature-coverage phase diagrams for various choices of Rt (from Ref. 0]) 
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FIG. 3: Groundstate phase diagrams for the Ising square lat- 
tice with pairwise and three-spin interactions: (a) i? = 1/4; 
(b) 7?= 1/2. 



phase transitions. Sample data for the 4th order cumu- 
lant and specific heat are shown in Fig. H] the specific 
heat peaks diverge with increasing lattice sizes for fields 
above H/\Jnn\ = —2.0, but for fields more negative and 
close to the c(2 x 2) phase boundary, the specific heat 
peaks first decrease and then diverge again with larger 
lattice sizes. Similar behavior was observed for the 4th 
order cumulant: for the small lattice, there is more than 
one diverging correlation length; but if the lattice sizes 
are big enough, only one dominates. Therefore, in certain 
range of small lattice sizes, the behavior is easy to con- 
fuse with XY-like 111, and GPU accelerated simulations 



of large lattice sizes is essential. For positive fields and 
low temperatures there is hysteresis in the m vs. H data 
indicating the presence of first-order transitions, but at 
higher temperature the data obtained for increasing and 
decreasing fields are essentially identical. 

The resultant phase diagram in field-termperature 
space is shown in Fig. [S] The c(2 x 2) phase is sepa- 
rated from the disordered phase on the high field side 
by a phase boundary which contains a tricritical point 
but on the low field side the transition appears to stay 
second-order down to the lowest temperature studied, 
and the row-shifted (2 x 2) state is also bounded by 
a line of second order transitions. Since the transi- 
tion from the c(2 x 2) ordered phase to the disordered 
phase should belong to the Ising universality, we expect 
ctjv — (logarithm) for a second-order phase transition, 
and a/i^ = 2 for a first-order phase transition. At the 
tricritical point, the exact(conjectured) value for the ex- 
ponent Oijv is |, which is supported by many renormal- 
ization group calculations [l2|. Therefore, we estimate 
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FIG. 4: Bulk properties for 7? = 1/4, Rt = -1/4: Specific 
heat c and 4tli order cumulant U of the corresponding order 
parameter for (a) H/\Jnn\ = 1- Data are for: L=32, ■; L=64, 
x; L=128, o; L=256, (b) H/\J„n\ = -5. Data are for 
L=30, ■; L=40, V; L=64, x; L=128, o; L=168, +; L=256, 



the exponent from finite size behavior of specific heat 
peaks, Cmax ~ L"/'^ , near the connecting section of the 
first- and second-order transition line, and found the tri- 
critical point is close to ksT /\Jnn\ — 0.592. To confirm 
our estimation and get more accurate location, we also 
calculate the density distribution of the order parameter, 
as shown in Fig. [5] The final estimation of the tricirtial 
point is A:bT/|J„„| = 0.5915(4), i7/|J„„| = 2.98073(8), 
and the evaluated exponent a/v = 1.59(2) agrees nicely 
with the predicted value. 

The corresponding phase diagram in coverage- 
temperature space is shown in Fig. [71 Here we see that 
the c(2 X 2) phase and the L.G.-|-c(2 x 2) coexistence 
phase, which is present below the tricritical point, ap- 
pear over substantial ranges of and T, whereas the 
row-shifted (2 x 2) phase is actually confined to a very 
narrow range of coverage. This phase diagram is sub- 
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FIG. 5: Phase diagram in magnetic field-temperature space 
for R = 1/4, i?i = —1/4. The sohd curves are second-order 
phase boundaries and the dashed line indicates first-order 
transitions. The open triangle indicates the location of a tri- 
critical point. 



120 r 




FIG. 6: Tricritical point(T=0.5915(4), H=2. 98073(8)) for 
R = 1/4, Rt — —1/4: (a) Curve fit of the specific heat peaks, 
(b) The density distribution of order parameter Mc(2x2) for 
L = 256. 
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FIG. 7: Temperature-coverage phase diagram for R 
1/4, Rt = -1/4. 
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FIG. 8: Adsorption isotherms of the lattice gas model with 
7? = 1/4, i?t = —1/4. The arrows mark the second-order 
phase transitions. 



stantially different from tliat predicted in Fig. and 
in particular there is no triple point. However, If third- 
nearest- neighbor two-body interactions are added, the 
ground state degeneracy for the (2 x 2) state will be re- 
moved. Then, tricritical points involving the (2 x 2) phase 
could occur, and triple points in the region of the first 
order transition perhaps as well, i.e., the predicted phase 
diagram in Fig. (^t could then be valid. 

In Fig. [8] we show adsorption isotherms which are 
obtained for several different temperatures and for 
comparison include the Langmuir isotherm which would 
be correct for a non-interacting lattice gas. The jump 
in the low temperature data shown by the dotted line 
clearly locates the first-order transition; but the second- 
order transitions, indicated by the arrows, are extremely 
difficult to identify from the adsorption isotherms. The 
step-like behavior of the lowest temperature adsorption 
isotherm shown in this figure is not dissimilar to the 
multiple risers which are seen for multilayer adsorption, 
but here it merely represents multiple transitions within 




FIG. 9: Field-temperature phase diagram for R = 1/2, Rt = 
— 1. The solid curves are second-order transitions and the 
dashed lines show first-order phase boundaries. The open 
square indicates the location of the terminating critical point. 



a single layer! 

B. R = 1/2, Rt ^ -1 

The same thermodynamic properties were determined 
as were described in Sec. A, and since there were no sig- 
nificant differences in the nature of the results, we shall 
not show any raw data for this case. The resultant phase 
diagram in H-T space is shown in Fig. ^ A line of first- 
order transitions, terminating in a critical point, sepa- 
rates a lattice liquid from a lattice gas state, and a line 
of second-order transitions bounds a row-shifted (2 x 2) 
phase. 

Since the lack of symmetry among the two different 
phases at the critical point that terminates the first order 
line, the relevant scaling fields r, h are comprised by lin- 
ear combinations of the thermodynamic fields T, H asfisj 



T = T-Tc + s{H - He 



h = H - Hc + r{T - Tc) 



(9) 



(10) 



where s and r are parameters controlling the extent of 
field mixing. As a result, the associated conjugate scaling 
operators £,A4 are also linear combinations of the spin- 
spin interaction energy density u and the magnetization 
m as 



£ = 
M = 



u — rm 

1 — rs 
m — su 



1 
1 

TV 



rs 



(11) 

(12) 
(13) 



= ^(^(JiCTj + R^aiak + Rt ^ cno-jcrk) (14) 
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FIG. 10: The density distributions for L = 64 and 128 at the 
critical point fcsr/|J„„| = 2.2738(4), ///1J„„ | = 1.24925(13). 
The full curve is the corresponding distribution for the two- 
dimensional Ising model for L = 400. 
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FIG. 11: Temperature-coverage phase diagram for 7? 
1/2, Rt = -1. 



where N is the total number of spins. According to the 
finite-size scaling [T3|. the joint probability distribution 
Pl{£, -M) near criticality should obey the following scal- 
ing ansatz: 

PL {£,M)- K+^K+pe,M (A+, JA^ , K+5£, AmK A^rXlS) 
where, 

As = a^L^I^ Km = aML"-^'" (16) 
A£A+ - A^A+ = (17) 
5M = M-<M>c, 6£ = £- <£ >c (18) 

The subscripts c denotes that the averages are taken at 
criticality. For appropriate choices of the nonuniversal 
factors as and qm, function ps^m would be universal. 
After integration over £, exactly at criticality, where h = 
r = 0, one has 



along the transition line, but the reduced exponents j/v 
and P/v seems to belong to Ising universality, similar to 
what we found in Ref Q. 

The phase diagram is replotted in coverage- 
temperature space in Fig. [TT] Here, too, we see 
that the row-shifted (2 x 2) phase is present only over 
a relatively narrow range of coverages, but the L.G. 
-|- L.L. coexistence phase is stable over a much larger 
region oi T — 9 space. Comparing with Fig. we see 
that there are qualitative differences between the actual 
behavior and the phase diagrams which had previously 
been "guessed"; but again, if the degeneracy allowing 
for row-shifted structures were removed, the predicted 
phase diagram in Fig. [2]i might hold. 



IV. CONCLUSIONS 



Pl{M) ajJ^L^/'^PUL^/'^ajJ^SM) (19) 

where the function p^ characterizes the universality 
class, the form of which has been well established for 
the two-dimensional Ising model. In Fig. [TOl we plot 
the density distribution function at estimated criticality 
kBT/\Jnn\ = 2.2738(4), iJ/|J„„| = 1.24925(13) with the 
controlling parameter s — —0.30(2) for i = 64 and 128. 
The superimposed curve is the corresponding distribu- 
tion for the two-dimensional Ising model for L = 400. 
The nonuniversal factors aj^ for each lattice sizes is cho- 
sen in such a way that the variable aJ^L^^'^{Ai— < 
Ai >c) has unit variance. 

As for the critical exponents for the continuous tran- 
sition from the row-shifted (2 x 2) phase to the param- 
agnetic phase, the correlation length exponent u changes 



Monte Carlo simulations have been used to extract 
phase diagrams for simple models on a square lattice 
with three-body interactions which are larger in magni- 
tude than those which have been previously studied. We 
find qualitatively different behavior than that which had 
been suggested in Ref Q . If third- nearest- neighbor two- 
body interactions are added, the ground state degeneracy 
for the (2x2) state will be removed; but of course there 
is no guarantee that the finite temperature behavior will 
not show remnants of this effect. These results further 
demonstrate the complexity which may be found in rel- 
atively simple models with competing interactions. This 
problem is of interest to statistical mechanics in its own 
right and we believe that further Monte Carlo studies of 
such models will continue to display many of the features 
observed in experimental studies of adsorbed monolayers. 
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